#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
#include <mpi.h>

#define TRUE    1;
#define FALSE   0;

typedef struct{
    int J_ORDER;
    int J_ROW_TEST;
    float J_ERROR;
    int J_ITE_MAX;
    float **MA;
    float *MB;
}Jacobi;

typedef struct{
    float *S_SOLUTION;
    float *S_SOLUTION_AUX;
    int S_ITE;
}Solucao;

void inicializa_Jacobi(Jacobi *J, FILE *file){
    int i, j;

    fscanf(file, "%d", &(J->J_ORDER));
    fscanf(file, "%d", &(J->J_ROW_TEST));
    fscanf(file, "%f", &(J->J_ERROR));
    fscanf(file, "%d", &(J->J_ITE_MAX));

    J->MA = (float**) malloc(J->J_ORDER * sizeof(float*));


    for(i=0; i < J->J_ORDER; i++)
        J->MA[i] = (float*) malloc(J->J_ORDER * sizeof(float));

    J->MB = (float*) malloc(J->J_ORDER * sizeof(float));

    for(i=0; i < J->J_ORDER; i++){
        for(j=0; j < J->J_ORDER; j++){
            fscanf(file, "%f", &(J->MA[i][j]));
        }
    }

    for(i=0; i < J->J_ORDER; i++)
        fscanf(file, "%f", &(J->MB[i]));
}

int converge(Jacobi *J){
    int i, j;
    int aux;

    for(i=0; i < J->J_ORDER; i++){
        aux = 0;

        for(j=0; j < J->J_ORDER; j++){
            if(i != j)
                aux += J->MA[i][j];
        }

        if(aux >= J->MA[i][i])
            return FALSE;
    }

    return TRUE;
}

void Jacobi_Richardson(Jacobi J, Solucao *S, int id, int num_proc, int tam_calc){
    int i, j, k;
    MPI_Status status;

    S->S_ITE = 0;
    S->S_SOLUTION = (float*) malloc(J.J_ORDER * sizeof(float));
    S->S_SOLUTION_AUX = (float*) malloc(J.J_ORDER * sizeof(float));




    if(id == 0){

	float rcv_erro1 = 0.0f;
	float rcv_erro2 = 0.0f;
	int iterac = 0;

	do{
		for(k = 1 ; k < num_proc ; k++){
			MPI_Bcast(&rcv_erro1,1,MPI_INT,k,MPI_COMM_WORLD);
			MPI_Bcast(&rcv_erro2,1,MPI_INT,k,MPI_COMM_WORLD);	
			MPI_Bcast(&iterac,1,MPI_INT,k,MPI_COMM_WORLD);
			S->S_ITE += iterac;
			printf("%d\n",S->S_ITE);
		}
	}while((rcv_erro1/rcv_erro2 > J.J_ERROR) && (S->S_ITE < J.J_ITE_MAX));

    }else{
			float aux1 = 0.0f, aux2 = 0.0f, erro1 = 0.0f, erro2 = 0.0f;
			int iteracao = 0;
			int ini;
			int fim;
			if(id == 1)
				ini = (id-1)*tam_calc;
			else
				ini = (id-1)*tam_calc+1;
				
			fim = (id*tam_calc);
	
			for(i=ini; i < fim; i++)
			    S->S_SOLUTION_AUX[i] = S->S_SOLUTION[i];

			omp_set_num_threads(3);
		
			#pragma omp parallel for private(j)
			for(i=ini; i < fim; i++){
			    S->S_SOLUTION[i] = J.MB[i];

			    for(j=ini; j < fim; j++){
				if(i != j)
				    S->S_SOLUTION[i] -= S->S_SOLUTION_AUX[j] * J.MA[i][j];
			    }

			    S->S_SOLUTION[i] /= J.MA[i][i];
			}

			erro1 = 0.0f;
			erro2 = 0.0f;
			for(i=ini; i < fim; i++){
			    aux1 = S->S_SOLUTION[i] - S->S_SOLUTION_AUX[i];
			    aux2 = S->S_SOLUTION[i];

			    if(aux1 < 0)
				aux1 *= -1;
			    if(aux2 < 0)
				aux2 *= -1;

			    if(aux1 > erro1)
				erro1 = aux1;
			    if(aux2 > erro2)
				erro2 = aux2;
			}
			iteracao = (int)S->S_ITE++;
			MPI_Send(&erro1,1,MPI_FLOAT,0,0,MPI_COMM_WORLD);
			MPI_Send(&erro2,1,MPI_FLOAT,0,0,MPI_COMM_WORLD);
			MPI_Send(&iteracao,1,MPI_INT,0,0,MPI_COMM_WORLD);

	}

}

void imprimi_Solucao(Jacobi J, Solucao S){
    int i;
    float result = 0.0f;


    for(i=0; i < J.J_ORDER; i++)
        result += S.S_SOLUTION[i] * J.MA[J.J_ROW_TEST][i];

    printf("Iterations : %d\nRowTest: %d => [%f] ?= %f\n", S.S_ITE, J.J_ROW_TEST, result, J.MB[J.J_ROW_TEST]);
}

int main(int argc, char *argv[]){
    Jacobi J;
    Solucao S;
    FILE *data;
    int id, num_proc,tam_calc;
    
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &id);
    MPI_Comm_size(MPI_COMM_WORLD, &num_proc);
    
    data = fopen("matriz250.txt","r");

    inicializa_Jacobi(&J, data);

    fclose(data);

    if(converge(&J)){
	
	tam_calc = J.J_ORDER/(num_proc-1);
        Jacobi_Richardson(J, &S, id, num_proc, tam_calc);
        imprimi_Solucao(J, S);
    }
    else
        printf("Sistema linear divergente!");

    MPI_Finalize();

    return 0;
}

